// DEFINITIONS
#define ERROR_RETURN(retval) { fprintf(stderr, "Error %d %s:line %d: \n", retval,__FILE__,__LINE__);    exit(retval); }

// INCLUDES
#include "matrix.h"
#include <papi.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>

// PROTOTYPES
void matrixMul(float*, float*, float*, int);

int
main(int argc, char** argv)
{
	printf("Running %s...\n", argv[0]);

	// Seed RNG.
	srand(time(NULL));

	// Process args.
	if(argc < 2)
	{
		printf("ERRR Not enough arguments.\n");
		return -1;
	}

	int width = atoi(argv[1]);

	if(width < 1)
	{
		printf("ERRR Need positive size argument: %d.\n", width);
		return -1;
	}

	if(width % 16 != 0)
	{
		printf("ERRR Need size to be multiple of 16: %d.\n", width);
		return -1;
	}

	// Initialize matrices.
	int size = width * width * sizeof(float);

	float* M = (float*) malloc(size);
	float* N = (float*) malloc(size);
	float* P = (float*) malloc(size);

	matrixRandInit(M, width);
	matrixRandInit(N, width);
	//matrixOnesInit(M, width);
	//matrixOnesInit(N, width);

	matrixZeroInit(P, width);

	// Set up PAPI.
	// Make sure PAPI library version is correct.
	char errstring[PAPI_MAX_STR_LEN];
	int retval;
	
	if((retval = PAPI_library_init(PAPI_VER_CURRENT)) != PAPI_VER_CURRENT )
	{
		fprintf(stderr, "Error: %s\n", errstring);
		return -1;
	}

	float rtime, ptime, mflops;
	long long flpins;

	// Start counting!
	//if( (retval=PAPI_start(EventSet)) != PAPI_OK)
	if( (retval=PAPI_flops(&rtime, &ptime, &flpins, &mflops)) != PAPI_OK)
		ERROR_RETURN(retval);

	// Run mmm.
	matrixMul(M, N, P, width);

	// Get counts.
	if( (retval=PAPI_flops(&rtime, &ptime, &flpins, &mflops)) != PAPI_OK)
		ERROR_RETURN(retval);

	// Close PAPI.
	PAPI_shutdown();

	// Print out results...
	printf("width:\t%d\nrtime:\t%f\nptime:\t%f\n", width, rtime, ptime);

	// Free malloc'd memory.
	free(P);
	free(N);
	free(M);
	
	printf("Done!\n");
	
	return 0;
}

void
matrixMul(float* M, float* N, float* P, int width)
{
	/*
	printf("width:\t%d\n", width);
	printMatrix(M, width);
	printMatrix(N, width);
	*/

	for(int i = 0; i < width; i++)
	{
		for(int j = 0; j < width; j++)
		{
			float sum = 0;
			
			for(int k = 0; k < width; k++)
			{
				float a = M[i*width + k];
				float b = N[k*width + j];
				sum += a*b;
			}

			P[i*width + j] = sum;
		}
	}

	//printMatrix(P, width);
}
